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Abstract 

Background: The individual character of pharmacokinetics is of great importance in the risk assessment of new 
drug leads in pharmacological research. Amongst others, it is severely influenced by the properties and inter- 
individual variability of the enzymes and transporters of the drug detoxification system of the liver. Predicting 
individual drug biotransformation capacity requires quantitative and detailed models. 

Results: In this contribution we present the de novo deterministic modeling of atorvastatin biotransformation 
based on comprehensive published knowledge on involved metabolic and transport pathways as well as 
physicochemical properties. The model was evaluated on primary human hepatocytes and parameter identifiability 
analysis was performed under multiple experimental constraints. Dynamic simulations of atorvastatin 
biotransformation considering the inter-individual variability of the two major involved enzymes CYP3A4 and 
UGT1A3 based on quantitative protein expression data in a large human liver bank (n = 150) highlighted the 
variability in the individual biotransformation profiles and therefore also points to the individuality of 
pharmacokinetics. 

Conclusions: A dynamic model for the biotransformation of atorvastatin has been developed using quantitative 
metabolite measurements in primary human hepatocytes. The model comprises kinetics for transport processes 
and metabolic enzymes as well as population liver expression data allowing us to assess the impact of inter- 
individual variability of concentrations of key proteins. Application of computational tools for parameter sensitivity 
analysis enabled us to considerably improve the validity of the model and to create a consistent framework for 
precise computer-aided simulations in toxicology. 



Background 

The discovery and development of new drug entities is 
strongly handicapped by the circumstance that about 
50% of the drug candidates fail in the clinical studies 
[1]. About one quarter of candidate drugs fail due to 
toxicity or pharmacokinetic (PK) problems [2], and cur- 
rently, it is a well known fact, that toxicity is the major 
cause of attrition in the drug development process [3]. 
Therefore, it is quite clear that dangerous properties of 
drug entities have to be revealed very early in the drug 
evaluation studies [4]. Despite the ever growing effort to 
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apply computational power towards improving the 
understanding and in-silico prediction of drug pharma- 
cokinetics and response, the overall impact on preclini- 
cal safety testing has been modest. 

Application of systems biology holds tremendous pro- 
mise because it aims to understand and quantitatively 
describe biological phenomena within the framework of 
the hierarchical levels of metabolic pathways and regula- 
tory networks at the different scales of cells, tissue, 
organs and ultimately whole organisms [5,6]. However, 
despite emerging consensus that such a holistic 
approach is essential to provide the framework of pre- 
dictive toxicology, the number of successful case studies 
is still minuscule [7-9]. 

Current activities can be grouped into 
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(1) quantitative structure-activity-relationship (QSAR) 
methods, computational models based on compound 
structure and focused on potential interactions of small 
molecules with major classes of proteins such as drug 
metabolizing enzymes [10-15], transporters [16] and 
receptors [16-18]. Also important are physicochemical 
properties of the drug, for example solubility and per- 
meability that are estimated from the molecular struc- 
ture [19-22]. 

(2) in vitro kinetics for prediction of in vivo drug 
clearance using kinetic data from recombinant cyto- 
chromes P450 (CYPs), microsomes and hepatocytes 
(IVIVE: in vitro-in vivo extrapolations) [23]. 

(3) physiologically based PK (PBPK) modeling [24-28] 
which considers the anatomical, physiological and che- 
mical aspects of ADME (absorption, distribution, meta- 
bolism and elimination of the drug) [29-31] in multi- 
compartment models [32]. 

In addition to these simulations based on mathemati- 
cal models various computational and bioinformatics 
approaches are applied to extract information from high 
throughput data of drug response experiments at cellu- 
lar, tissue, organ and whole organism level. 

A critical assessment of the aforementioned tools, 
essentially to outline gaps that must be addressed for 
more reliable predictive simulation-based toxicology, 
indicates needs for more rigorous network models 
focusing at systems dynamics beyond kinetics of indivi- 
dual enzymes, consideration of inter-individual variabil- 
ity and systematic investigations of parameter sensitivity 
and its impact on model verification, discrimination and 
reduction, to name a few. The first issue is related to 
the design of the dynamic models for the drug elimina- 
tion process in the hepatocyte, which should be based 
on the integration of membrane transport processes for 
substrates and products as well as phase I and phase II 
reactions. These models need to be tightly linked to sti- 
mulus (dose) -response experiments. 

The issue of model parameterization in the context of 
modeling in toxicology has been already addressed in 
1995 by Andersen et al [24]. In addition to the problem 
of identifiability, particular attention should be given to 
correlation between parameters, very common in biolo- 
gical systems. 

Yet another question of interest concerns the subtle 
integration of the enormous inter-subject variability in 
enzymatic phenotypes into the model. This is of outer- 
most importance for predictions in toxicology and also 
in clinical pharmacology in order to design optimal 
treatments for individual patients. Consideration of this 
variability in phenotype should rest on quantitative pro- 
teomics and activity data from human liver tissue or 
hepatocytes representing a statistically significant por- 
tion of the population. 



The importance of this issue has been extensively 
addressed in the review of Rostami-Hodjegan and 
Tucker [23] and discussed in the context of IVIVE (in 
vitro-in vivo extrapolation) approach. Significant pro- 
gress has been made by this group through a strategy 
which is based on the simulation of drug disposition in 
virtual liver populations. The basic idea of this method 
is to portray the variability into the r max /K M -va\ues 
estimated from in vitro studies using human liver 
microsomes, primary or cryopreserved hepatocytes or 
recombinant expressed enzymes. This is a promising 
approach; however, improvement of the strategy is ima- 
ginable by separation of the two parameters r max and 
K M (maximal enzymatic rate and enzymatic affinity 
parameter in Michaelis-Menten- and derived kinetics). 

This concept, which is a focus of the present article, is 
driven by the possibility to incorporate separate infor- 
mation from pharmacokinetics and quantitative proteo- 
mics as well as future incorporation of the regulatory 
network responsible for variation of expression level of 
the enzymes, transporters and receptors. As model drug 
we chose atorvastatin (AS), one of the most frequently 
used 3-hydroxymethylglutaryl coenzyme A reductase 
inhibitors (statins). Statins thereby reduce cholesterol 
synthesis and they also stimulate the uptake of LDL- 
cholesterol from the blood. Although they are relatively 
safe drugs, the lipid-lowering effect of statins is inade- 
quate in some patients, and unpredictable drug-drug 
interactions can occur, as well as hepatic and extrahepa- 
tic adverse effects including hepatotoxicity, myopathy, 
and rare but severe rhabdomyolysis. 

Some aspects of dynamic modeling of statins have 
been tackled in previous studies, including the determi- 
nistic modeling of transport kinetics [33,34] and the 
pharmacokinetic modeling of the clearance mechanisms 
including consideration of unspecific binding effects 
[35]. In the latter case, a multi-compartmental model 
approach has been derived on isolated rat hepatocytes. 
However, the metabolism of AS inside the liver cell has 
not been described in details, missing the kinetic 
description of the different metabolite formations, which 
thus far precludes the consideration of inter-individual 
variability, apart from the fact that rat hepatocytes had 
been investigated. 

This study describes a deterministic modeling 
approach of the dynamic biotransformation and trans- 
port processes of AS in human hepatocytes considering 
both kinetics of transport processes as well as intracellu- 
lar detoxification processes via Phase I and Phase II 
enzymes. Attention is given to the detailed enzymatic 
and chemical reactions, including conversions between 
the acid and lactone form of AS as well as unspecific 
binding between metabolites and macromolecules like 
proteins. The additional integration of quantitative 
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protein data of Phase I and Phase II enzymes enables 
the dynamic analysis of the underlying inter-individual- 
variability of expression levels of these enzymes. 

Methods 

Chemicals 

Williams Medium E (WME) without phenol-red and 
without L-glutamine and stabilized L-glutamine were 
obtained from Pan-Biotech-GmbH (Aidenbach, Ger- 
many). Penicillin/Streptomycin and ITS-X were obtained 
from Invitrogen (Karlsruhe, Germany). Bovine serum 
albumin (BSA) and dimethyl sulfoxide (DMSO) were 
purchased from Sigma-Aldrich Chemicals (Taufkirchen, 
Germany). Dexa Inject was obtained from Jenapharm 
(Jena, Germany). AS, its metabolites and deuterated 
standards were purchased from Toronto Research Che- 
micals Inc. (North York, Canada). PBS was obtained 
from Invitrogen (Karlsruhe, Germany) and Complete 
Mini, EDTA-free from Roche Diagnostics (Mannheim, 
Germany). NaPP, sucrose and LiChrosolv were pur- 
chased from Merck (Darmstadt, Germany); acetonitrile 
from Roth (Karlsruhe, Germany), formic acid from 
Fluka, (Germany), UGT1A3 monoclonal mouse antibody 
from Abeam (Cambridge, England). Trypsin was pur- 
chased from Promega (Mannheim, Germany). Synthetic 
peptides were purchased from Sigma-Genosys (Haver- 
hill, UK). 

Isolation and cultivation of primary human hepatocytes 

Tissue samples from human liver resections were 
obtained from patients undergoing partial hepatectomy. 
Experimental procedures were performed according to 
the guidelines of the charitable state-controlled founda- 
tion HTCR (Human Tissue and Cell Research) Regens- 
burg, Germany, and the institutional guidelines for liver 
resections of tumor patients with primary or secondary 
liver tumors, Technical University Munich, MRI, 
Munich, Germany. The use of human hepatocytes for 
research purposes was approved by the local ethics com- 
mittees of the Ludwig-Maximilians-University of 
Munich [36] and the Charite, Humboldt University Ber- 
lin [37], Germany, and written informed consent was 
obtained from all patients. Hepatocytes were cultured 
on collagen gel precoated 6-well plates at a density of 
1.5-10 6 cells/well. Cells were allowed to attach to the 
collagen layer. After transport, culture media was dis- 
posed and attached cells were cultured 24 h at 37°C in a 
humidified chamber with 95%/5% air/C02 in serum-free 
medium WME, supplemented with albumin (0.1% (v/v)), 
penicillin/streptomycin (100 U/ml), stabilized L-gluta- 
mine (2 mM), dexamethasone dihydrogenphosphate 
(0.025% (v/v)) and ITS-X (5 mg insulin, 3.35 \ig 
natrium-selenit, 2.75 mg transferrin and 1 mg ethanola- 
mine), further named SFM. 



Time-series experiments 

Incubation with AS was started by disposing the culture 
media and cultivation of the attached cells at 37°C in a 
humidified chamber with 95%/5% air/C0 2 in 2 ml SFM, 
supplemented with 10 [iM AS, 0.1% (v/v) BSA and 0.1% 
DMSO. At specified time-points, three wells were 
further treated for the preparation of samples for the 
measurement of extracellular and intracellular metabo- 
lites, respectively. SFM media was collected and 50 [iL 
formic acid and deuterated internal standard was added 
for the further measurement of extracellular metabolites. 
Cells were harvested in pre-cooled albumin-free SFM, 
disrupted by freeze/thaw and ultra-sonification and cen- 
trifuged. The supernatant was used for the determina- 
tion of intracellular metabolites. 

For the preparation of samples for the protein mea- 
surements, culture medium was disposed and cells were 
harvested in pre-cooled PBS, supplemented with Com- 
plete Mini EDTA-free (1 Tablet/10 ml Buffer). Cell sus- 
pensions were centrifuged 5 min (500 g) at 4°C and cell 
pellets were resuspended in 150 \iL NaPP-buffer (0.1 M, 
pH 7.4), containing 250 mM sucrose and Complete 
Mini EDTA-free (1 Tablet/10 ml Buffer). Cells were dis- 
rupted by ultra-sonification and lyophilized for the ana- 
lysis of total protein concentration, CYP3A4 and 
UGT1A3 content. 

For cell number determination, cells from two wells 
were fixed with methanol-acetic acid fixative solution 
(10 min at 37°C and 4°C) and afterwards nuclei were 
stained for 15 min with Meyers Hamalaun (Sigma- 
Aldrich Chemie GmbH, Germany), rinsed with water 
and air-dried. Stained nuclei were counted in digital 
images (10 per well) at 40-fold Magnification (ImageJ 
Image Processing and Analysis Program). 

Quantification of atorvastatin and its metabolites 

AS and ASL, and their para- (ASpOH, ASLpOH) and 
ortho-hydroxy-metabolites (ASoOH, ASLoOH), were 
determined by LC-MS-MS analysis using the respective 
deuterium labeled analogues as internal standards, 
essentially as described [38]. HPLC separation was per- 
formed at 30°C on a XBridge Shield RP18 column (2.1 
x 50 mm, 3.5 |im, Waters) using (A) 1 mM formic acid 
and (B) acetonitrile as mobile phases at a flow rate of 
0.4 ml/min. Gradients were programmed as follows: 
63% A for 4 min; linear decrease of A to 60% within 9 
min; linear decrease of A to 55% within 2.5 min; 55% A 
for 1 min; increase of A to 63% in 0.2 min. Equilibration 
time of the column was 20 min. MS-MS analysis was 
performed on an Esquire HCT ultra ion trap mass spec- 
trometer (Bruker Daltonics, Bremen, Germany) coupled 
to an HPLC 1100-System (Agilent, Waldbronn, Ger- 
many) consisting of binary pump G1312A, degasser 
G1379A, well-plate sampler G1367A and column 
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thermostat G1330B. The ionization mode was electro- 
spray (ESI), polarity positive, mass range mode ultra- 
scan, and nitrogen was used as a drying and nebulizer 
gas. The following parameters were applied: nebulizer 
45 psi, dry gas 10 1/min, dry temperature 300°C, capil- 
lary 4100 V, scan range 200 - 600 m/z. 

Precursor and product ions (m/z) of analytes and inter- 
nal standards, respectively, were ATV (559 and 440.2; 
466.2), [ 2 H 5 ]ATV (564 and 445.2; 471.2), ATV-L (541.2 
and 448.2), [ 2 H 5 ]ATV-L (546.2 and 453.2), p-OH-ATV 
(575 and 440.2; 466.2), [ 2 H 5 ]p-OH-ATV (580 and 445.2; 
471.2), p-OH- ATV-L (557 and 448.4), [ 2 H 5 ]p-OH-ATV-L 
(562 and 448.4), o-OH-ATV (575 and 466.4), [ 2 H 5 ]o-OH- 
ATV (580 and 471.2), o-OH-ATV-L (557 and 448.4), 
[ 2 H 5 ]o-OH- ATV-L (562 and 448.4). Sample quantifica- 
tion was possible in a range from 0.5 to 500 pmol. 

CYP3A4 and UGT1A3 protein quantification 

Protein quantification of CYP3A4 and UGT1A3 in 
human liver microsomes and relative protein quantifica- 
tion of UGT1A3 in lyophilized samples of primary hepa- 
tocytes was performed by immunoblotting as described 
previously [38,39]. 

Cell lysates for absolute quantification analysis of 
CYP3A4 were prepared from lyophilized human primary 
hepatocytes by sonification in the presence of glass 
beads in buffer containing Complete Mini-Protease Inhi- 
bitor Cocktail, and following homogenization. In the 
CYP3A4 quantification assay, three synthetic isotopically 
labeled peptides (13C/15N amino acid) were used as 
internal standard for calibration. These isotope labeled 
standard peptides represented sequence analogues to 
proteotypic peptides of CYP3A4, which arise from tryp- 
tic digestion. After acetone precipitation and resolving 
of the proteins in 8 M Urea, a definite amount of inter- 
nal standard peptides was added. The sample mixture 
was reduced with 5 mM DTT and alkylated with 15 
mM iodacetamide in 50 mM ammonium bicarbonate. 
Subsequently samples were digested with trypsin at 42°C 
for 4 h (enzyme/substrate ratio of 1:10). Efficiency of 
tryptic digestion was checked by SDS-PAGE followed by 
silver staining. The resulting peptides were purified 
using CI 8 OMIX® Tips (Varian, Darmstadt, Germany) 
according to manufacturer's suggested protocol and 
separated on a nanoliter-flow Ultimate HPLC system 
(Dionex, Idstein, Germany). After injection (15 |il), pep- 
tides were trapped and desalted on a precolumn (0.3 
mm I.D. x 5 mm PepMapTM, Dionex) at a flow rate of 
30 (il/min in 0.1% TFA for 6 min. Peptides were trans- 
ferred to the separation column (75 [im I.D. x 250 mm 
PepMapTM column, Dionex) and separated in a linear 
gradient of mobile phase (A: 0.1% formic acid, B: 84% 
acetonitrile/0.1% formic acid) from 5% B to 35% B over 
a period of 35 min with a flow of 290 nl/min. The 



column effluent was continuously directed into the 
NanoSpray II source of a 4000QTrap mass spectrometer 
(Applied Biosystems, Foster City, CA, USA). The MS 
was set up to run a multiple reaction monitoring experi- 
ment essentially, as described previously [40], including 
two to three parent-to-product ion transitions for each 
internal standard peptide as well as the corresponding 
transitions of native peptide of CYP3A4. The instrument 
settings were as follows: ion spray voltage, 3-4 kV; inter- 
face heater, 150°C; declustering potential, 50 V; collision 
energy, peptide specific; entrance potential, 10 V; colli- 
sion cell exit potential, 10 V. MS data were processed 
by integrating the appropriate peak areas from extracted 
ion chromatograms by MultiQuantTM Software 
(Applied Biosystems). The absolute amount of CYP3A4 
protein was calculated from the peak area ratio "internal 
standard peptide/native peptide". Total protein content 
of the samples were determined by amino acid analysis 
(AAA) on a Waters 2695 HPLC system using the 
AccQ»Tag derivatization method (Waters, Eschborn, 
Germany), according to manufacturer's instructions. 

Identification of metabolic network structure 

Numerous metabolic and physicochemical aspects about 
AS had to be considered in the initial model building. 
AS exists in two forms, a very lipophilic lactone (ASL) 
and a comparably hydrophilic hydroxyl-acid (AS). AS is 
converted enzymatically via an instable intermediate 
product into ASL, mediated by different UGT isoen- 
zymes [41,42]. Recent investigations by ourselves and 
others have shown that the most important contributor 
to UGT-driven lactonization is UGT1A3, whereas 
UGT1A1 plays an insignificant role [38,42]. AS acids 
and lactones are inter-converted chemically into each 
other [43]. However, studies have indicated, that the 
chemical lactonization of AS to ASL can be neglected at 
physiological pH of 7.4 [43]. Recent studies highlight, 
that different PON enzymes might also be possible con- 
tributors to the lactone hydrolysis and that PON1 is 
present in liver [44-48]. 

Both AS and ASL are hydroxylated in human hepato- 
cytes leading to para- and ortho-hydroxy-metabolites, 
ASpOH, ASoOH, ASLpOH and ASLoOH [49,50], 
mainly catalyzed by CYP3A4 [51]. Recent studies have 
reported, that CYP2C8 and CYP3A5 also hydroxylate 
AS to a minor extent [50,52]. 

Furthermore, AS is transported into the cell via 
organic anion transport polypeptides (OATP), and 
recent studies on recombinant systems showed, that 
OATP1B1 and OATP2B1 contribute to the AS import 
[53,54], which are both expressed in the human liver 
[55,56] 

OATP1B3 is supposed to be also a main contributor 
to drug transport [54], since it shows high gene 
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expression levels in liver [57], but its importance for AS 
has not been investigated kinetically so far. 

OATP transporters have been reported to be bidirec- 
tional facilitated diffusion transporters, independent 
from ATP and Na+, K+ and H+, but with a possible 
involvement of reduced glutathione [58,59]. However, 
previous investigations on transport mechanisms on rat 
hepatocytes, showed, that the intracellular concentra- 
tions of pitavastatin and other compounds are much 
higher than outside the cell [33,34,60]. Facilitated or 
passive diffusion do not allow a greater intracellular 
than extracellular concentration of the parent drug of 
interest, when the source is the initial extracellular 
concentration, because it is a concentration-gradient 
dependent mechanism. Therefore, the import mechan- 
ism should be rather considered as active mediated 
transport, which is not concentration-gradient depen- 
dent, rather than as the proposed facilitated diffusion 
process. 

As shown with the recombinantly expressed transpor- 
ters mentioned above [54], transport of the acidic meta- 
bolites, ASpOH and ASoOH, by OATP1B1 was similar 
to that of AS, and transport of the corresponding lac- 
tones was also mediated by this transporter, although at 
somewhat lower rates. We therefore assumed that the 
same OATP transporters are responsible for the import 
of AS and its hydroxylated and lactone metabolites into 
hepatocytes. 

AS and its hydroxylated metabolites, ASpOH and 
ASoOH, are actively exported out of the hepatocytes 
into the bile by the ATP-dependent MDR1 transporter 
[61,62]. In addition, acidic and lactone form of AS 
showed inhibitory effects in transport studies of sub- 
strates of MDR1 and MRP2 [63-65], pointing to the 
competitive transport mechanism of these substrates at 
this proteins. Further, the transporters MRP1, MRP3 
and MRP6 are also reported to be responsible for the 
transport of organic compounds including AS from 
inside the hepatocytes into the plasma [56,66]. 

Passive diffusion might play also an important role in 
the transfer of AS, ASL and the corresponding metabo- 
lites, ASpOH and ASoOH, ASLpOH and ASLoOH, 
respectively. Since the acidic forms of AS are rather 
hydrophilic and the lactone forms of AS are rather lipo- 
philic, it can be assumed that passive diffusion plays a 
more important role for the lactone forms, and the 
transporter mediated active transport plays a more 
important role for the acidic forms, as reported earlier 
for statins [67]. 

Finally, lipophilic drugs have a high affinity to bind 
non-specifically to proteins, and previous studies con- 
centrated on the modeling of drug binding in the intra- 
cellular and the extracellular space as well as on the 
surface of the cells [60,68,69]. 



Mathematical modeling 

The mathematical model can be described by the system 
of non-linear ordinary differential equations 



d (cjVcomp) 



dt 



(1) 



where the change in extracellular, intracellular or 
unspecific bound metabolite concentration c ; in the 
extra - or intracellular compartment with V comp is 
effected by the conversion or production of contributing 
chemical or enzymatic reactions or by transport steps 
Tip respectively. The extracellular volume, VZmp> equals 
to the volume of the media used. The intracellular 
volume, Vcomp> equals to the total volume of all cells 
used, and is determined by multiplying the cell number 
by the volume of a single hepatocyte, estimated to be 
14.1 pL by the approximation of a spherical shape with 
a diameter of 30 [im [70]. 

Appropriate reaction kinetics r, ; are modeled for the 
CYP3A4 hydroxylation, the UGT1A3 lactonization, the 
chemical and enzymatic lactone hydrolysis and intracel- 
lular unspecific binding to macromolecules. 

Previous studies determined substrate inhibition 
kinetics of the CYP3A4 mediated hydroxylation of AS 
on human microsomes [52]. However, inhibition effects 
contribute severely only at a concentration higher than 
100 (iM. Furthermore, our model approach considers 
the competitive nature of alternative substrates, by inte- 
grating the CYP3A4 hydroxylation of AS and ASL as 
reaction kinetics describing the competitive conversion 
of alternative substrates to alternative products, 
illustrated for the hydroxylation of AS to ASpOH 



7max,3A4, ASpOH 



CAS 



^3A4,ASpOH = 



KM,3A4ASpOH 



Den ■■ 



CAS 



1 + Den 

CAS 



(2) 



KM,3A4,ASpOH Km,3A4,ASoOH 

casl casl 



Km ,3 A4 ,ASLpOH Km,3A4,ASLoOH 

(see Additional file 1 - Derivation of Atorvastatin 
kinetics at CYP3A4). 

The lactonization of AS to ASL is mediated by 
UGT1A3 enzymes and the reaction is formulated as 
substrate inhibition kinetics [41]. 



TlA3,AS 



:,1A3,AS ' C AS 



Km,1A3,AS + C AS + 



(CAST 



(3) 



I,1A3,AS 



The lactone metabolites are either hydrolyzed chemi- 
cally to the respective acid metabolites [43] inside (c) or 
outside (m) the cell, or enzymatically by the contribution 
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of PON enzymes inside the cell. Both reactions are 
described as first order kinetics 



7 c/m -i c 

Thydj = "CR ' j + KPON,j ' Cj 



(4) 



Unspecific binding of intracellular metabolites to 
macromolecules is formulated as 



»*wrf-**((^- 1 ) < j- < j)' 



(5) 



with the dissociation coefficient kjis and the intracellu- 
lar fraction unbound [71] 



fUj 



c c 



C j,eq + C j,eq 



(6) 



which describes the ratio between the intracellular free 
concentration c c - to the sum of intracellular free and 
bound concentration, c c - and in equilibrium (index eq) 
(see Additional file 2 - Derivation of kinetics of unspeci- 
fic protein binding). However, the intracellular free and 
bound concentrations in equilibrium are not measur- 
able; therefore, the fraction unbound fuj is set as a para- 
meter to be estimated in the parameter optimization 
procedure. 

Transport steps include active import and export of 
the metabolites as well as passive diffusion steps. Both 
active import by OATP1B1 or OATP2B1 and export 
of AS are described as Michaelis-Menten-kinetics 
[53,54] 



Tim/ ex, AS 



7max,lBl/2Bl/ex,AS * CAS 
KM,lBl/2Bl/ex,AS + C AS 



(7) 



whereas the active transport kinetics of the other 
metabolites are assumed to be of first-order [72]. 



im/ex,) ~ kim/ex,) ' Cj 



c/m 



(8) 



Besides the active transport, metabolites undergo pas- 
sive diffusion through the double-layer lipid-membrane. 
Passive diffusion, described as 



Jiff 



D, 



j-A cells -(cf-<f), 



(9) 



is driven by the concentration difference between 
outside and inside the cell, (c™ — (f\ over the lipid- 
membrane with thickness d, through all cells with the 
total surface area A ceUs , and controlled by the diffusion 
coefficient Dj and is comprised as the permeability 
coefficient J?u 



Optimization procedure for estimation of model 
parameters 

The optimization procedure is based on evolutionary 
strategies which are implemented with JavaEva (WSI 
Computer Science Department, Center for Bioinfor- 
matics, University of Tubingen, Germany) and a MVA 
(main vector adaptation) mutation operator [73]. The 
optimization procedure estimates the parameters in 
equations (2) to (9) based on the optimization criterion 



arg mm 

P 



N } 

EEi 

n=l j=l \ 



I r calc 



(10) 



where the deviation of calculated and measured con- 
centrations divided by the measurement standard devia- 
tion Sp squared and summed over all metabolites J and 
all time points N, has to reach a minimum. Additional 
optimization constraints 2Lrefu AS >fu AS u Pasl >Pas> Pas 
>Pasoh and Pasl >Pasloh> because the lactones have a 
higher lipophilicity than the acids and the metabolites 
are supposed to be more hydrophilic than the corre- 
sponding parent lactone or acid drug. 

The integration of the differential equations (1) using 
the reaction kinetics in equations (2) to (9) was performed 
by the differential algebraic equation solver LIMEX (Kon- 
rad-Zuse-Zentrum fur Informationstechnik, Berlin) [74] . 

Relative abundance approach for prediction of r max - 
parameters 

For pharmacokinetic predictions taking inter-individual 
variability of CYP3A4 and UGT1A3 expression levels 
into account, maximal rate parameters are predicted via 
a relative abundance approach, which is based on the 
assumption that the maximal rate of the reaction is pro- 
portional to the enzyme concentration: 



max 



max, reference, j 



4 



reference 



(ii) 



The maximal velocities r max ,i,j of the respective 
enzyme e, here CYP3A4 or UGT1A3, in the conversion 
to the product Pj in the liver of interest li is estimated 
from the respective maximal rate r m ax, reference,; an d the 
enzyme concentration ^ e f erence in the reference liver and 
the enzyme concentration cf in the liver of interest li. 

Computational approach 

The mathematical model of AS metabolism was coded 
in FORTRAN language and linked to the numerical 
integrator LIMEX, also written in FORTRAN. After 
compilation to the executable program, optimization 
was started by the call of JavaEva. The mathematical 
model is supplemented as SBML-file for review purpose. 
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Results 

Model-setup 

The model structure of AS biotransformation in primary 
human hepatocytes is schematized in Figure 1 and com- 
prises the description of the extracellular, intracellular 
and unspecifically bound metabolites, as well as corre- 
sponding reaction and transport steps and unspecific 
protein binding. 

Since the calibration and quality samples were prepared 
in the same media as the experimental metabolite samples 
and because it can be assumed that the intracellular pro- 
tein concentration is much higher compared to the out- 
side of the cells, the unspecific protein binding was only 



considered for the intracellular space. Regarding the unspe- 
cific protein binding, it can be assumed that the intracellu- 
lar protein concentration is much higher compared to the 
outside of the cells. This assumption is based on the esti- 
mation of the ratio of intracellular protein to extracellular 
media protein concentration. From the total cell protein 
measurement and the determination of total cell volume a 
total intracellular protein concentration of about 30 g/1 
can be estimated. The media includes as sole protein com- 
ponent 0.1% (v/v) albumin, corresponding 1 g/1. Therefore, 
this assumption is justified. 

Because our hepatocyte culture model does not allow 
to distinguish experimentally import and export, single 





OATP1B1 
OATP2B1 




Unspecific 
intracellular binding 
to macromolecules 





MDR1 




MRP2 





Figure 1 Scheme of the full model of Atorvastatin metabolism in primary human hepatocytes. The scheme illustrates atorvastatin 
metabolic and transport pathways in human hepatocytes. The model comprises the active import and export of AS and ASL, and corresponding 
para- and ortho-hydroxy-metabolites, ASpOH and ASoOH, ASLpOH and ASLoOH, mediated by import proteins, OATP1B1 and 0ATP2B1 (solid 
lines), and export proteins, MDR1 and MRP2 (solid lines), as well as passive diffusion steps (dashed lines). AS and ASL are hydroxylated to the 
corresponding metabolites by CYP3A4, respectively (solid lines). Compound AS is converted via an instable glucuronid-intermediate to ASL 
mediated by UGT1A3 (solid lines). The lactones are hydrolyzed either chemically or enzymatically by putative PON to the respective acids (solid 
and dotted line). Acids and lactones are further subject to unspecific binding to macromolecules in human hepatocytes (dashed lines). 
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irreversible, actively mediated transport steps for both 
directions, as well as passive diffusion steps were consid- 
ered for each compound, respectively, similar to 
mechanistic transport models of other studies [33,34]. 
Only in the case of AS different transport steps for 
OATP1B1 and OATP2B1 were implemented in the 
model, because specific KMs had been determined in 
previous studies [53,54]. 

Model verification 

Model verification was performed by quantitative mea- 
surements of AS and all metabolites in the extra- and 
intracellular space during time-series experiments, per- 
formed first on a single batch of primary hepatocytes 
from a single individual 1. Special emphasis was given 
to the issue of parameter sensitivity and correlation 
between parameters to improve the quality of the 
model. 

Time series/stimulus-response measurements on primary 
human hepatocytes 

From the metabolite concentrations in primary human 
hepatocytes of individual 1 (see Additional file 3 - Ator- 
vastatin metabolite concentrations from the time-series 
experiment on primary human hepatocytes of individual 
1) we calculated total recovery at each time point (Table 
1). The recovery is calculated from material balance 
equations and is defined as the sum of intracellular and 
extracellular metabolite amounts at the respective time- 
point divided by initial AS amount. As evident from 
Table 1, the total recovery was close to 100% after ten 
minutes but decreased at higher time-points. This may 
be explained by unspecific protein binding leading to 
the pool of bound metabolites, which increases over 



Table 1 Recovery of atorvastatin metabolites in the time- 
series experiment on primary human hepatocytes of 
Individual 1 



time [min] 


Recovery [%] 


0 


100 


10 


97.6 


30 


96.8 


60 


94.8 


120 


88.8 


180 


86.5 


240 


89.0 


300 


89.0 


360 


80.8 


480 


n.o. 


600 


78.2 



The recovery, calculated from material balance equations at each 
experimental time-point, is defined as the sum of intracellular and 
extracellular metabolite amounts divided by measured initial atorvastatin acid 
amount (n.o.: not observed) (see Additional file 3 - Atorvastatin metabolite 
concentrations from the time-series experiment on primary human 
hepatocytes of individual 1). 



time due to intracellular accumulation of metabolites. 
This result points to the necessity and importance of 
the implementation of unspecifically bound metabolites 
in the model 
Parameter estimation 

Parameter optimizations were performed with aid of 
evolutionary algorithms (JavaEva, [i = 8 parent and X = 
4 children) and using the nominal parameter values of 
the reactions, transport and diffusion steps in equations 
(2) to (9). The optimization criterion in equation (10) 
was used, which demands to minimize the difference 
between experimental concentration data and model 
simulation. In the optimization procedure certain para- 
meters were fixed to values (Table 2), which were either 
identified in or assumed from previous investigations. 
The Michaelis-Menten constants K M of the CYP3A4 
hydroxylation and of the AS transporters OATP1B1 and 
OATP2B1 were fixed to values determined with recom- 
binant enzymes [50,53,54]. The K M and Kj of the 
UGTlA3-lactonization were fixed to values determined 
on human liver microsomes [41], The rate constant of 
spontaneous hydrolysis of ASL, I<cr> was estimated from 
experimental observations [43] and assumed to be same 
for the hydrolysis of the lactone metabolites, ASLpOH 
and ASLoOH, respectively. The dissociation rate con- 
stant kdis of unspecific binding was fixed to a value, con- 
sidered to be very high in a previous study on modeling 
of protein binding mechanisms [33]. 
Analysis of the predicted concentration-time-profiles 
The model predicted concentration-time-profiles are 
illustrated together with the measured concentrations 
(see Additional file 3 - Atorvastatin metabolite concen- 
trations from the time-series experiment on primary 
human hepatocytes of individual 1) in Figure 2. Both 
intracellular AS and ASL are converted to the corre- 
sponding para- and ortho-hydroxy metabolites, ASpOH 
and ASoOH, and ASLpOH and ASLoOH, respectively. 
However, the acidic metabolites, ASpOH and ASoOH, 

Table 2 Parameters which were fixed in the optimization 
procedure to literature values 



Parameter Value Units Literature 



K|Vi,3A4,ASpOH 


25600 


pmol-ml" 1 


[50] 


Km,3A4,ASoOH 


29700 


pmol-ml" 1 


[50] 


KM,3A4,ASLpOH 


1400 


pmol-ml" 1 


[50] 


Km,3A4,ASLoOH 


3900 


pmol-ml" 1 


[50] 


Km,1B1,AS 


18900 


pmol-ml" 1 


[54] 


K)V1,2B1,AS 


200 


pmol-ml" 1 


[53] 


Km,1A3,AS 


12000 


pmol-ml" 1 


[41] 


K|,1A3,AS 


75000 


pmol-ml" 1 


[41] 


k C R 


0.0025 


min" 1 


[43] 


kdis 


600 


min" 1 


[33] 



Parameter values were adopted from literature sources outlined in the right 
column. 
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Figure 2 Measured concentrations and simulated fits of Atorvastatin compounds, individual 1. Displayed are the measured 
concentrations (rhombs) with mean and standard deviation (n = 3) of Atorvastatin compounds in the culture medium (upper two rows) and 
free Atorvastatin compounds inside primary hepatocytes (lower two rows) of individual 1, and corresponding simulation profiles (solid lines). First 
and third row: AS (left), ASpOH (middle) and ASoOH (right). Second and fourth row: ASL (left), ASLpOH (middle) and ASLoOH (right). 



show higher intracellular concentrations than the 
lactone metabolites, ASLpOH and ASLoOH. The ratio 
between intracellular AUC 0 -6oo min of the acidic form to 
that of the lactone form equals 20.1 in case of the para- 
hydroxy-metabolites and 23.6 in the case of the ortho- 
hydroxy-metabolites. 

Accordingly, the extracellular concentrations of the 
acidic metabolites are higher than that of the lactone 
metabolites. The ratio between extracellular AUC 0 , 600 
min of the acidic form to that of the lactone form equals 
54.3 in case of the para-hydroxy- metabolites and 70.8 in 
the case of the ortho-hydroxy-metabolites. 

Further, the results show that the components have 
higher concentrations in the intracellular space than 
outside the cell. The ratio of intracellular AUC 0 _ 600 min 
to extracellular AUC 0 _ 600 min equals minimally 4.5 in the 
case of ASoOH and maximally 38.2 in the case of ASL- 
pOH. Similarly, the ratio of intracellular c max to extra- 
cellular c max equals minimally 2.5 in the case of ASoOH 
and maximally 34.5 in the case of ASLpOH. 



Simultaneous model verification on different individual 
hepatocyte donors 

The process of model verification so far has been 
applied on a single experiment on primary human hepa- 
tocytes of individual 1. In the next step it has to be pro- 
ven, that the model is also capable to describe different 
individual metabolic profiles, especially being further 
able to reflect the inter-individual parameter variability, 
not only in the parameters r max of the phase I and II 
reactions, but also in the transporters. Therefore, the 
model verification is performed based on AS biotrans- 
formation data on primary hepatocytes from three dif- 
ferent individuals simultaneously. 

Therefore, maximal rate constants r max of CYP3A4 
hydroxylation and UGT1A3 lactonization were predicted 
for individual 2 and 3 via the relative abundance 
approach (equation (11)) using the r m ^_parameters of 
individual 1 (Table 3) and the protein concentrations 
observed on the primary human hepatocytes (Table 4). 
These parameters were then fixed in the optimization. 
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Table 3 Verified parameters of CYP3A4 and UGT1A3 in 



individual 1 


Parameter 


Value 


rel. Error [%] 


Units 




r max,3A4,ASpOH 


1108 


5.5 


pmol-min" 1 


•ml" 1 


r max,3A4,ASoOH 


3345 


2.3 


pmol-min" 1 


•ml" 1 


r max,3A4,ASLpOH 


1228 


12.7 


pmol-min" 1 


■ml" 1 


r max,3A4,ASLoOH 


2756 


11.6 


pmol-min" 1 


•ml" 1 


r max,1A3,AS 


956 


14.8 


pmol-min" 1 


•ml" 1 



Parameters are listed with nominal values from optimization and relative 
parameter errors from FIM based identifiability analysis (see Additional file 4 - 
Analysis of parameter sensitivity and following model reduction procedure) in 
the corresponding units shown. 



Further, the rate constants k of PON mediated lactone 
hydrolysis of individual 2 and 3 were estimated in the 
optimization procedure. 

However, the analysis of parameter sensitivity and 
identifiability (see Additional file 4 - Analysis of para- 
meter sensitivity and following model reduction proce- 
dure) showed that the first order kinetics of OATP1B1 
mediated import is favored over the zero-order kinetics 
of OATP2B1 import of AS, and that the first order 
kinetics are hard to distinguish from first order passive 
diffusion mechanisms in this evaluation system. Thus, in 
the following, mechanisms of active transport and pas- 
sive diffusion are lumped, resulting in the apparent 
transport rate for import, 

Yim,j = kimj • Cj + Pj • Cj = K{m,j • Cj , (12) 

and in the apparent transport rate for export 

Tex,j = kex,j ' Cj + Pj ' Cj = ^ex,j ' Cj , (13) 

described by the product of apparent rate constant 
Kim/ex of import or export and extra- or intracellular 
concentration cj 1 ^, respectively (Figure 3). Consequently, 
the rate constants n of import and export are set indivi- 
dually, and were to be estimated in the optimization 
procedure. 

The simultaneous model verification was performed 
based on the stimulus response data obtained from pri- 
mary human hepatocytes of individual 1 (see Additional 
file 3 - Atorvastatin metabolite concentrations from the 

Table 4 Protein concentrations of CYP3A4 and UGT1A3 



in primary human hepatocytes 



Individual 


CYP3A4 [pmol ml' 


'] 


UGT1A3 [-] 


1 


1027 ±107(n 


= 2) 


1.00 


2 


611 ±120(n 


= 4) 


0.29 


3 


755 ±19(n 


= 2) 


0.10 



Protein concentrations were calculated from total CYP3A4 and UGT1A3 
content per sample divided by total intracellular volume of hepatocytes, 
respectively. Protein concentrations of CYP3A4 are listed in mean and 
standard deviation in the units pmol per ml intracellular volume. Relative 
protein concentrations of UGT1A3 are normalized to individual 1. 



time-series experiment on primary human hepatocytes 
of individual 1), individual 2 (see Additional file 5 - 
Atorvastatin metabolite concentrations from the time- 
series experiment on primary human hepatocytes of 
individual 2) and individual 3 (see Additional file 6 - 
Atorvastatin metabolite concentrations from the time- 
series experiment on primary human hepatocytes of 
individual 3), respectively. 

The model prediction is in satisfying agreement with 
experimental data of individual 2 (Figure 4). However, 
in case of individual 3, the model prediction is compar- 
ably poor, because there are relatively high deviations, 
especially with the intracellular metabolites ASL and 
ASoOH and extracellular ASpOH (Figure 5, solid line). 
One reason could be the extended contribution of beta- 
oxidation in AS metabolism. Beta oxidation at the 
heptanoic acid side chain is a typical transformation 
pathway for all statins, but is reported to play only a 
minor role in humans [75]. However, high activity of 
beta-oxidation of fatty acids, which is responsible for the 
supply of ATP for gluconeogenesis in type 2 diabetes 
mellitus [76,77], which individual 3 was diagnosed with, 
may severely influence the AS metabolism. To test this 
hypothesis, respective reactions with the acid metabo- 
lites as substrates were considered in the model veri- 
fication of individual 3, but the results showed no 
significant improvement (data not shown). The second 
reason could be, that CYP3A4 and UGT1A3 protein 
concentrations of individual 3 used in the estimation 
of corresponding r max value via relative abundance 
approach (equation (11)) differ from the measured mean 
value (Table 4). Due to the fact that the r max parameters 
of CYP3A4 and UGT1A3 show high sensitivities in the 
parameter sensitivity analysis (see Additional file 4 - 
Analysis of parameter sensitivity and following model 
reduction procedure), a variation in the values would 
have a high impact on the metabolic profiles. Therefore, 
in a further optimization step, the CYP3A4 protein con- 
centrations are allowed to vary in the interval of mean ± 
standard deviation (Table 4), where the standard devia- 
tion of UGT1A3 protein concentration is assumed to be 
30% of the mean value. Notably, an improvement in the 
model prediction on individual 3 could be achieved in 
the case of the intracellular metabolites AS and ASL 
and extracellular AS (Figure 5, dashed lines). But there 
are still major deviations in case of intracellular ASL 
and ASoOH and extracellular ASpOH, which could not 
be explained any further, but shows, that there must be 
some other reactions or effects in the system, which 
have not been discovered yet. 

The estimated model parameters for individual 1, 2 
and 3, summarized in Table 5, display the proposed 
parameter variability of enzyme mediated reactions 
as well as of the transport steps. However, the rate 
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Figure 3 Scheme of the model used in the simultaneous model verification. In contrast to the full model illustrated in Figure 1, this model 
contains lumped transport kinetics, which means, one import and one export rate for each metabolite, respectively. This transport simplification 
resulted from the local parameter sensitivity analysis (see Additional file 4 - Analysis of parameter sensitivity and following model reduction 
procedure). 



constants of ASLoOH import and export of individual 2 
and of ASL import and export of individual 3 are con- 
siderably higher than the other transport rate constants, 
which points to strong linear dependencies of these 
parameters. 

The mathematical model of AS metabolism in human 
hepatocytes of individual 1 is supplemented as SBML- 
file (see Additional file 7 - Model of atorvastatin meta- 
bolism in primary human hepatocytes of individual 1, 
and BioModels database). 

Dynamic analysis of inter-individual CYP3A4 and UGT1A3 
expression level variability 

Based on the model version of optimized parameters, 
gained in the simultaneous model fit, the effect of inter- 
individual variability of CYP3A4 and UGT1A3 protein 



expression levels was investigated by linking the protein 
expression data of 150 liver samples (Figure 6) via the 
described relative abundance approach in equation (11), 
using individual 1 as reference. UGT1A3 and CYP3A4 
protein concentrations of individual 1 were converted 
from based on total protein amount to based on micro- 
somal protein amount by multiplying with the factor 
0.22, determined in human liver homogenates and cor- 
responding microsome fractions via Bradford test [78]. 

Simulations were performed with an initial extracellu- 
lar AS concentration of 50 (pmol ml" 1 ), which is in the 
range of physiological plasma concentrations [54,79], 
over a time period of 1200 min. 

The most important question that arises from this 
dynamic analysis was, how the metabolic profiles of the 
intracellular metabolites AS, ASpOH and ASoOH are 
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Figure 4 Measured concentrations and simulated fits of Atorvastatin compounds, individual 2. Displayed are the measured 
concentrations (rhombs) with mean and standard deviation (n = 3) of Atorvastatin compounds in the culture medium (upper two rows) and 
free Atorvastatin compounds inside primary hepatocytes (lower two rows) of individual 2, and corresponding simulation profiles (solid lines). First 
and third row: AS (left), ASpOH (middle) and ASoOH (right). Second and fourth row: ASL (left), ASLpOH (middle) and ASLoOH (right). 



influenced by this variability, since they are considered 
to be the active drugs, which inhibit HMGCoA-reduc- 
tase [80]. Therefore, AUC, c max and t(c max ) of the con- 
centration-time-profiles of either AS alone or the sum 
of concentration of AS, ASpOH and ASoOH were cal- 
culated for each liver sample over a time period of 1200 
min and the distributions over all liver samples were 
evaluated, respectively. Finally, appropriate probability 
density functions are fitted to the distributions (Figure 7). 
The probability density function characteristics are sum- 
marized in Table 6. Obviously, there are differences 
between the examination of only AS or the sum of con- 
centrations of AS, ASpOH and ASoOH. AUC, c max and 
t(c max ) have lower values in case of AS alone compared 
to the sum of all acidic metabolites. The population 
mean of AUC is 75051 (pmol ml" 1 min) in the case of 
AS and 203617 (pmol ml" 1 min) in the case of the sum 
of AS, ASpOH and ASoOH. Also, c max is lower in the 
case of AS alone, 201 (pmol ml" 1 ), compared to the sum 



of the acidic metabolites, 366 (pmol ml" 1 ). Further, the 
maximal concentration appears at a shorter time point, 
48 min, in the case of AS alone, compared to the time 
point, 100 min, of the sum of AS, ASpOH and ASoOH. 
The results are quite explainable, since ASpOH and 
ASoOH are the hydroxylated products of AS and there- 
fore their maximal concentrations event at a delayed 
time-point compared to AS. 

Further, except for the AUC of AS, the probability den- 
sity functions show a very narrow shape, regarding the 
relative standard deviation, when comparing to the sample 
standard deviations of the underlying distribution of 
CYP3A4 and UGT1A3, respectively. The probability den- 
sity functions fitted to CYP3A4 and UGT1A3 distributions 
have a relative standard deviation of 259% and 137%, 
respectively, whereas the relative standard deviations of 
the probability density functions of AUC, c max and t(c max ) 
are lower than 50%, except for the probability density 
function of AUC of intracellular AS, which equals 123%. 
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Figure 5 Measured concentrations and simulated fits of Atorvastatin compounds, individual 3. Displayed are the measured 
concentrations (rhombs) with mean and standard deviation (n = 3) of Atorvastatin compounds in the culture medium (upper two rows) and 
free Atorvastatin compounds inside primary hepatocytes (lower two rows) of individual 3, and corresponding simulation profiles (solid lines and 
dashed lines). First and third row: AS (left), ASpOH (middle) and ASoOH (right). Second and fourth row: ASL (left), ASLpOH (middle) and ASLoOH 
(right). The dashed lines show the additional consideration of putative beta-oxidation of atorvastatin acids, AS, ASpOH and ASoOH, and of 
adaptation of CYP3A4 and UGT1A3 protein concentrations in the model fit. 



Discussion 

The deterministic modeling of drug metabolism has a 
major advantage compared to traditional pharmacoki- 
netic models. The distinction between metabolism and 
elimination inside the hepatocytes and the detailed 
description of the biotransformation network structure 
allows the observation of the reaction kinetics and the 
transport mechanisms involved. Our modeling consid- 
ered the major acid and lactone metabolites in the 
extracellular and intracellular space as well as appropri- 
ate reaction kinetics and transport steps. An experimen- 
tal limitation in primary hepatocytes concerns the lower 
limit of quantification of each metabolite. Therefore, a 
relatively high initial concentration of AS of about 10 
uM had to be chosen in order to permit measurement 
of the described major metabolites. Interestingly, the 
intracellular concentration of the parent drug AS was 
higher than the extracellular concentration. This 



indicates, that either the OATP mediated import does 
rather follow an active uptake mechanism than the pro- 
posed facilitated diffusion mechanism [58,59], or so far 
unknown transporters may be involved in the transport 
steps. By any means it justifies the chosen implementa- 
tion of kinetics for the active mediated transport beside 
the passive diffusion mechanism. In contrast to concen- 
trations determined in plasma in clinical studies 
[54,79,81], extracellular concentrations of the acidic 
metabolites were much higher compared to the lactones. 
A possible explanation of these discrepancies may be 
that the in-vitro investigation on isolated primary 
human hepatocytes was performed over rather short 
time intervals (hours), whereas analyses in humans 
usually cover longer periods (days). A second explana- 
tion could be the relatively high initial concentration of 
AS of about 10 uM used in this study, compared to 
plasma concentrations observed to be lower than 200 
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Table 5 Parameter Variability in the simultaneous model 
verification on individual atorvastatin metabolism data. 



Parameter 


Individual 


Individual 


Individual 


1 Inifr 

Units 




1 


2 


3 






1 108 


660 


606 


pmol'iTiin ^ -ml ^ 


ASpOH 










^max,3A4, 


3345 


1991 


1830 


pmol-min" 1 -mr 1 


ASoOH 










^max.SA^ 


1228 


731 


672 


pmol-min" 1 -mr 1 


ASLpOH 










r max,3A4, 


Z/DO 


1 04U 


I jUo 


-1 -1 
pmol-min -ml 


ASLoOH 










fm ~, v 1 AC 

max, i f\D,f\D 


957 


281 


120 


pmol-min" 1 -mr 1 


' "KUIM.AjL 


308 


99.2 


0.00 


10" 3 min -1 


^PON ASLOH 


280 


545 


82.9 


10" 3 min" 1 


^ex,AS 


2.1 7 


0.16 


0.1 7 


1 1 1 -min ^ 


^ex,ASL 


21.8 


2.48 


10.5 


uL-min ^ 


^ex,ASpOH, 


0.80 


1.25 


0.52 


uL-min 1 




1.60 


3.62 


0.91 


uL-min' 1 


^ex,ASLpOH 


1.13 


1.01 


0.59 


uL-min 1 


^ex.ASLoOH 


2.67 


86.3 


0.77 


uL-min -1 


^im.AS 


20.3 


4.42 


3.45 


uL-min -1 


^im.ASL 


275 


25.8 


82.7 


uL-min -1 


^im,ASpOH 


3.96 


2.01 


2.53 


uL-min -1 


^im.ASoOH 


0.39 


0.22 


0.00 


uL-min -1 


^ex.ASLpOH 


33.7 


7.24 


8.59 


uL-min -1 


^im.ASLoOH 


26.1 


530 


2.43 


uL-min -1 


^B-OxAS 






15.0 


10" 3 min" 1 


l<B-Ox,ASOH 






0.00 


10" 3 min" 1 



listed for the individuals 1, 2 and 3 with nominal values from 
the corresponding units shown. 



Parameters are 
optimization in 



(pmol ml" 1 ). Furthermore, we found that the recovery of 
metabolites decreases over the time. This confirms on 
the one hand the contribution of unspecific binding to 
macromolecules, most severely in the intracellular space, 
as observed previously [35], but on the other hand 
could also be contributed to a certain extent by the 
unspecific binding to the collagen layer or the plates 
[34], especially of highly lipophilic ASL and lactone 
metabolites. However, due to the washing procedure 
preliminary to the cell harvesting and disruption proce- 
dure, this effect could not be observed in the chosen 
experimental set-up. 

The parameter identification was considered satisfac- 
tory, as the simulation profiles fitted well to the corre- 
sponding measured metabolite concentration data. But, 
it was questionable if the optimized parameters could be 
considered to be sensitive and identifiable. Therefore, a 
local parameter sensitivity and identifiability analysis 
based on the Fisher-Information-Matrix was performed. 
The results showed that the full model of AS metabo- 
lism is not identifiable. Consequently a model reduction 
procedure was set-up and could be applied successfully, 
indicating that the parameter identification difficulties 
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Figure 6 Individual variability of CYP3A4 and UGT1A3 protein 
concentration level. Displayed is the distribution of CYP3A4 (top) 
and UGT1A3 (bottom) protein concentration data of individual 
human liver microsomes (n = 150). Protein concentrations are 
normalized to minimal value, respectively. 



are caused by non-identifiable transporter steps, dis- 
played by low parameter sensitivities and high parameter 
errors and correlations, which were then reduced in the 
model. However, this outcome does not necessarily 
mean that the non-identifiable transport steps are not 
present or not used in human hepatocytes, but rather 
implies, that they cannot be distinguished from the 
remaining transport steps in the model verification. 
Thus, the remaining transport steps in the model cap- 
ture the probable superimposed contribution of several 
transport mechanisms in-vivo. 

Unfortunately, the local identifiability of the final 
model version reveals some uncertainties caused by 
remaining high correlations present in the transport 
steps, as well as in the intracellular reaction network. 
An additional problem attributed to the modeling and 
parameter optimization is related to the fact that several 
parameters have lumped characteristics, because they 
are used for more than one compound, like the fraction 
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Figure 7 Distributions from dynamic analysis of the Inter- 
individual variability in Atorvastatin metabolism. Distributions 
(bars) of AUC (top), c max (middle) and of t(c max ) (bottom), and fitted 
probability density functions of simulated profiles of AS (solid lines) 
and of the sum of AS, ASpOH and ASoOH (dashed lines), 
respectively. Distributions are predicted by the dynamic simulation 
of Atorvastatin metabolism, implementing individual CYP3A4 and 
UGT1A3 protein concentration data (Figure 6). 



unbound factor, which is set to be the same for all 
acidic or lactone AS compounds, respectively. A possible 
solution would be to set individual parameters for each 
compound. However, then the parameter identifiability 
difficulties would be even more severe. This is due to 
the fact, that the experimental observation on primary 
human hepatocytes is constrained strongly on several 
limitations. For example, the range of initial concentra- 
tions is strongly limited to the quantification of metabo- 
lites of interest. Further, the number of data points is 
strongly limited to the available cell number of primary 
hepatocytes, isolated from surgery removal. Finally, also 
the individual character of human cells due to individual 



Table 6 Probability density functions fitted to inter- 



individual variability distributions 



Probability density 
function 


Mean 


s.d. 


rel. s.d. 
[%] 


AUC [pmol ml ^min] 










AS 


Log-logistic 


75051 


92038 


123 


AS+ASpOH+ASoOH 


Log-logistic 


203617 


35775 


18 


Cmax [pmol ml" 1 ] 










AS 


Log-logistic 


201 


64 


32 


AS+ASpOH+ASoOH 


Log-logistic 


366 


10 


3 


t(c max ) [min] 










AS 


Log-logistic 


48 


18 


39 


AS+ASpOH+ASoOH 


Log-logistic 


100 


12 


12 



Probability density functions fitted to distributions of AUC, c max and t(c max ) 
(depicted in Figure 7), generated in dynamic simulations considering inter- 
individual variability of CYP3A4 and UGT1A3 protein concentration data 
(Figure 6). Distributions of AUC, c max and t(c max ) were fitted best with Log- 
logistic function. Displayed are the mean, the standard deviation (s.d.), and 
the relative standard deviation (rel. s.d.) of the fitted function, respectively. 

genetic and medical background precludes an exact 
reproduction of experiments and the reproducibility of 
experimental results. The remaining parameter uncer- 
tainty seems to be crucial for the prediction of the 
population probability on the first sight. Actually, this is 
not the case. Since the parameter errors and the correla- 
tions indicate which parameter changes in what extent 
do not influence the computed time courses, the varia- 
bility in the respective parameter error range would not 
influence the probability distributions. In general, the 
question if parameter errors and insensitivities are suffi- 
ciently small can only be answered when knowing the 
requirements of the pharmacological application as 
drug, which was not given in this study. 

The simultaneous model verification based on primary 
human hepatocytes from different individuals illustrates 
the individual character of drug metabolism in the liver. 
Further, it indicates the inter-individual variability of rate 
parameters not only of the phase I and II enzymes, but 
also of the transport enzymes. In case of patients with 
type 2 diabetes, AS metabolism is probably influenced 
strongly by beta-oxidation and further so far unknown 
effects. Therefore, the individual model verification might 
also be used to reveal undesired pharmacokinetic effects. 
However, this analysis tool should be checked carefully 
on other drug systems in the future, too. 

The individual nature of the AS metabolism was investi- 
gated in this study in the domain of inter-individual varia- 
bility of CYP3A4 and UGT1A3 protein expression levels 
in human hepatocytes, by performing dynamic analysis on 
the verified model, individualized by implementation of a 
comprehensive set of protein data from a liver bank. The 
results show, that the inter- individual variability strongly 
affects the biotransformation behavior and therefore 
reveals the individual character of AS metabolism. 
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Consequently, this individuality in the pharmacokinetics of 
AS also points out the individual character of pharmaco- 
dynamics at the drug target, namely the HMGCoA-reduo 
tase in human hepatocytes. 

However, this subject- variability can also be expected 
to be present in the transport protein expression, as 
previously shown for OATP-C [82]. Therefore, this 
variability should be taken into account in the future by 
implementing corresponding population expression data 
in the dynamic analysis. Another source of variation 
is the well-known polymorphism in CYP- and UGT- 
enzymes and likely transport proteins, which very prob- 
ably causes inter-individual variations of the catalytic 
activity and the substrate affinity K M . Thus, further 
efforts should also consider this individual difference, 
leading to the improvement of the model prediction. 

To resolve the contributing transport mechanisms in 
more detail, additional information should be imple- 
mented from transport investigations on recombinant 
systems, as described previously [72,83], which enable 
the differentiated identification of both basolateral and 
apical transporters of AS and the corresponding trans- 
port kinetics. Further, also the estimation of unspecific 
drug protein binding could be improved in the future by 
using radiolabeled compounds in the experiments and 
modeling approaches as described previously [60,69]. 

Conclusions 

We believe that the results of our simulations provide 
strong arguments for rigorous dynamic modeling of drug 
biotransformation at the cellular level embedded in a sys- 
tems biology approach. By resolving the detailed metabolic 
network structure with metabolites and catalyzing 
enzymes, we investigated the dynamic variation of atorvas- 
tatin metabolism affected by the inter-individual variability 
of expression levels of phase I and phase II enzymes. 

In contrast to experimental investigations on recombi- 
nant systems or tissue fractions of hepatocytes, like 
microsomes, the investigation on primary human hepa- 
tocytes enabled the holistic and most realistic in-vitro 
observation of drug biotransformation, because it is pos- 
sible to observe the coupled contribution of metabolism 
and transport to the entire processes. De novo of this 
study, we identified intracellular concentration profiles 
of atorvastatin metabolites in primary human hepato- 
cytes in a time-series approach. 

Such an approach is essential for integration of 
further-reaching issues, such as drug-drug interactions, 
impact of regulation networks linked to nuclear recep- 
tors and particularly to quantitatively account for sub- 
ject-variability. The integration of this variability caused 
by genetic or environmental variations is crucial for pre- 
dictive pharmacokinetic modeling. 



Such a rigorous modeling approach critically depends 
upon a tight link between experimental observations 
and model design, simulation and verification. While the 
results are promising, some limitations in the parameter 
identification were still encountered. In the long term, 
these open problems can only be solved by stronger 
links to other research areas, such as pharmacogenetics, 
characterization of transporters, etc. On the whole, our 
contention is that the problem of parameter identifiabil- 
ity is an indispensable ingredient of model verification. 
Systematic investigations - if possible linked to optimal 
experimental design - can greatly strengthen the cred- 
ibility of the models. 

However, we present a model that goes much further. 
The domain of application does not remain the system 
behavior for which it has been elaborated. The model 
provides the possibility to link further modules such as 
gene regulation, drug target metabolism and present the 
important links to be implemented into the PBPK envir- 
onment. Finally, the model structure used in this study 
should be considered as a module to be integrated into 
the framework of multi-scale whole body modeling and 
simulations necessary to tackle the drug disposition in 
patient populations. 

Additional material 



Additional file 1: Derivation of atorvastatin kinetics at CYP3A4. This 
file contains the derivation of the atorvastatin kinetics at CYP3A4, which 
describes the competition of alternative substrate degradation to 
alternative products (supplemented as .pdf-file). 

Additional file 2: Derivation of kinetics of unspecific binding to 
macromolecules. This file contains the derivation of kinetics of the 
unspecific binding of atorvastatin metabolites to macromolecules, for 
example proteins (supplemented as .pdf-file). 

Additional file 3: Atorvastatin metabolite concentrations from the 
time-series experiment on primary human hepatocytes of 
individual 1. Extracellular concentrations (upper part) and intracellular 
concentrations (lower part) of atorvastatin acid and lactone (AS and ASL) 
and corresponding para- and ortho-hydroxy-metabolites (acids: ASpOH 
and ASoOH; lactones: ASLpOH and ASLoOH) at the defined time-points 
with mean and standard deviation (n = 3) from triplicate measurements 
per LC-MS/MS (n.d.: not determinable; n.o.: not observed). The recovery is 
calculated from material balance equations and is defined as the sum of 
intracellular and extracellular metabolite amounts at the respective time- 
point divided by initial AS amount (supplemented as .pdf-file). 

Additional file 4: Analysis of parameter sensitivity and following 
model reduction procedure. This file contains the analysis of parameter 
sensitivity and the exemplary model reduction procedure, which is 
necessary for achieving a high quality predictive model of atorvastatin 
metabolism (supplemented as .pdf-file). 

Additional file 5: Atorvastatin metabolite concentrations from the 
time-series experiment on primary human hepatocytes of 
individual 2. Extracellular concentrations (upper part) and intracellular 
concentrations (lower part) of atorvastatin acid and lactone (AS and ASL) 
and corresponding para- and ortho-hydroxy-metabolites (acids: ASpOH 
and ASoOH; lactones: ASLpOH and ASLoOH) at the defined time-points 
with mean and standard deviation (n = 3) from triplicate measurements 
per LC-MS/MS (n.d.: not determinable) (supplemented as .pdf-file). 



Bucher et al. BMC Systems Biology 201 1, 5:66 
http://www.biomedcentral.eom/1752-0509/5/66 



Page 1 7 of 1 9 



Additional file 6: Atorvastatin metabolite concentrations from the 
time-series experiment on primary human hepatocytes of 
individual 3. Extracellular concentrations (upper part) and intracellular 
concentrations (lower part) of atorvastatin acid and lactone (AS and ASL) 
and corresponding para- and ortho-hydroxy-metabolites (acids: ASpOH 
and ASoOH; lactones: ASLpOH and ASLoOH) at the defined time-points 
with mean and standard deviation (n = 3) from triplicate measurements 
per LC-MS/MS (n.d.: not determinable) (supplemented as .pdf-file). 

Additional file 7: Model of atorvastatin metabolism in primary 
human hepatocytes of individual 1. The sbml-file contains the model 
of Atorvastatin metabolism in human hepatocytes with the model 
parameters identified on primary human hepatocytes of individual 1 
(supplemented as .xml-file). 



List of abbreviations 

AS: atorvastatin acid; ASL: atorvastatin lactone; ASOH: hydroxy-atorvastatin 
acids (para- and ortho-); ASpOH: pa ra-hydroxy-atorva statin acid; ASoOH: 
ortho-hydroxy-atorvastatin acid; ASLOH: hydroxy-atorvastatin lactones (para- 
and ortho-); ASLpOH: para-hydroxy-atorvastatin lactone; ASLoOH: ortho- 
hydroxy-atorvastatin lactone; B-OX: beta-oxidation of atorvastatin acids; c: 
index: intracellular (cytosol); calc: index: calculated; CR: index: chemical 
reaction of hydrolysis of lactones; CYP: cytochrome P450 monooxygenase; 
Den: Denominator; e: index: enzyme; eq: Index: equilibrium; ex: index: export; 
FIM: Fisher-Information-Matrix; fu: fraction unbound; /': index: reaction; im: 
index: import; / J; index: compound; k: index: parameter; rate constant; K M : 
Michael is-Menten constant; /: index: parameter; m: index: extracellular 
(medium); meas: index: measured; n:N; index: time-point; MDR: multidrug 
resistance protein; MRP: multidrug resistance-related protein; OATP: organic 
anion transport protein; PON: paraoxanase; r: reaction rate; rmax: maximal 
reaction rate; UGT: UDP-glucuronosyl-transferase; 1A3: index: UGT1A3; 1B1: 
index: OATP1B1; 2B1: index: OATP2B1; 3A4: index: CYP3A4; 
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